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AN ALTERNATING-DIRECTION-IMPLICIT ALGORITHM 
FOR THE UNSTEADY POTENTIAL EQUATION 
IN CONSERVATION FORM 

Richard R. Chiproan 
Grumman Aerospace Corporation 

1 - SUMMARY 

An implicit finite-difference scheme is presented for the efficient com- 
putation of unsteady potential flow about airfoils. The formulation uses density 
and velocity potential as dependent variables, and is cast in conservation form 
to assure the theoretically correct determination of shockwave location and 
speed. To enable boundary conditions to be imposed directly on the airfoil 
surface, a time-varying sheared-rectllinear coordinate transformation is em- 
ployed. Calculated time- history solutions on a pulsating airfoil are compared 
with the results of another unsteady transonic code. The method is demon- 
strated to ha\e excellent numerical stability and to give accurate solutions with 
sharply resolved shocks. 


2 - INTRODUCTION 

The transonic flow regime has long been known to be the most critical for 
flutter and other unsteady aeroelastic phenomena. Until recently, there was 
no efficient method for calculating unsteady aerodynamics in this speed range; 
consequently, transonic flutter prediction has relied on wind tunnel testing. 
With the advent of faster computers and the emphasis on transonic cruise and 
maneuver capabilities for new aircraft design, much progress has been recently 
made in the development of both steady and unsteady transonic computational 
methods. 

In unsteady transonic aerodynamics, work has proceeded along two dis- 
tinct lines. In the first, researchers have produced linear'zed unsteady solu- 
tions about nonlinear mean (steady) flows. The efforts of Ehlers^; Traci, 
Albano, and Farr*; Cunningham’; Liu^; and Fung, Yu, and Seebass’are 
examples of this approach. From experimental measurements, such as those 
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of Tijdemen*, it has been obvious that these llnearixed solutions are only valid 
for a limited set of problems. Consequently, other researchers have pursued 
a second approach - the use of finite-difference methods to obtain solutions to 
the coupled steady /unsteady flow. In this area, the works of Ma^us and 
Yoshihara^; Lerat and Sides'; Beam and Warming*; Ballhaus and Steger^'; 
Ballhaus and Goor]ian^^~^^; Isogai^'; Chipman and Jameson^*; Gooijian^*; 
and Sankar and Tassa^' are notable. The first three dted efforts in coupled 
steady /unsteady flow have produced methods for solving the full Euler equa- 
tions, which (although computationally too expensive for routine use) do pro- 
vide excellent benchmark calculations. Ballhaus' works have produced an 
efficient method for solvi.ig the lew -frequency, small-perturbation form of the 
potential equation , thus making possible economic solutions to a range of im- 
portant transonic unsteady problems. 

To extend this range, the latter four works have solved the full-potential 
flow equations. Isogai developed the first such procedure. The next two 
works are improvements in that: (1) conservation-law form is used to accurate- 

ly locate shock waves in space and time; and (2) the ADI scheme is used for 
computational efficiency. The Chipman -Jameson method uses density and ve- 
locity components as dependent variablis, giving a simple system of first-order 
equations to be solved. The Goorjian method uses velocity potential as the 
dependent variable and, by time-linearizing the density, derives a single 
(though complicated) scalar equation. The method of Sankar and Tassa uses 
the strongly-implicit-procedure algorithm and has been formulated both in con- 
servative and nonoonservative form; thus far, the scheme has been coded only 
in nonconservative form. 

The present method uses both the density and velocity potential as de- 
pendent variables, resulting in a simple system of two equations. It is a signifi- 
cant improvement over the authors' previous method in that the use of the 
potential function totally eliminates numerically created vorticity present in the 
prior method. Furthermore, it retains the desirable features of conservation 
form and efficient implicit differencing. 

The theoretical development of the present method is given in Section 4 of 
this report. The discussion covers the basic flow equation, the details of the 
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altern«tlnG'>Qireclton-iinplicit (ADI) algorithm, the addition of artifidal-viscotity 
terms to capture shocks, the introduction of a time-varying sheared- rectilinear 
coordinate transfor.nation , the specification of the boundary conditions, and the 
Implementation of the method for a spedflc test problem - an airfoil pulsating 
in thickness. 

In Section 5, the test problem is discussed in detail, and computed resxilts 
are presented. Thesp results are compared with those of Reference 15. Also 
studies of the effect of variations in the time-step size and artifldal viscosity 
are summarized. 

Section 6 presents conclusions and recommendations. 

3 - LIST OF SYMBOLS 


D.( ) 

D,( ) 

{f(T,{w})} 

h 

J 

M 

( )“ 
rl, r2 
S(x,t) 
t, T 


B local speed of sound 

> central difference in x of enclosed quantity 
^ central difference In y of enclosed quantity 
■ vector function of T and the vector W 

> stagnation enthalpy 

B Jacobian of transformation 
= Mach nun.ber 

B n-th time level of enclosed quantity 
= right-hand side of Eq (9) or Eq (26) 

= instantaneous airfoil-surface location 
- time 


u, U » streamwise component of velocity, 4 -aligned 

contravariant velocity component 

V, V * stream-normal component of velocity, rj-allgned 

contravariant velocity comi>onent 

{W} B general vector (in this case, having components p and d>) 

x,y B untransformed coordinates 

2 (x,y,t) * surface of flow discontinuity 

y B specific heat ratio 
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A( ) ■ time Increment of enclosed quantity 

•. ■ artificial vlscoBlty coefficient 

(,T) ■ transformed coordinates 

M ■ see Eq (£) 

p ■ density 

^ B velocity potential 


^ ^m 
( )- 
A(‘) 

B partial derivative of enclosed quantity with respect to "m" 
B free-stream value of enclosed quantity 
B solution to Eq (11) 

variable 


4 - THEORY 



Basic Flow Equations 


For unsteady, two-diitionslonal, isentropic potential flow, the equations of 
conservation of mass and momentum can be written 


Pf + (p4>,),-‘-(P0,), * 0 

(1) 


+ hB 0, 

(2) 

where 




(3) 


(Equation (2) is called the inte^ated unsteady Bernoulli equation.) Both 
equations are in conservation form, albeit (2) is considered weak-conservation 
form. Hence, (1) conserves mass across shocks, where the flow variables are 
discontinuous, and results in the jump condition 

[plrt + (pu]z, + [pv)Zy«0, (4) 

where z(x,y,t) = 0 represents the surface of discontinuity and [ ] denotes 
the jump of the enclosed variables across the discontinuity. 

Implicit Algorithm 

To simplify the ensuing discussion, (1) and (2) are represented by the 
vector equation 
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( 5 ) 


A family of difference schemes for solvin^f this system is 

{aw} + m At {f ■♦'} ♦ ( 1 - M ) At{f •} - 0 , (6) 

Os/i SI, 

where n and n 1 denote the present and new time levels. For u = 0, the 
scheme is explicit; otherwise it is implicit. The case p = 1/2 is the standard 
Crank-Nicolson scheme, which has second-order accuracy. Using a Taylor 
series expansion, the variables at the new time can be written 


{f>}«{f“}-K^{f}At + 0(At^). 

«{f“} + V-{AW}{f"}+0(At*). (7) 

Substitution into (6) gives 

{aw} + p Atr V - {aw} ){£“}=- At{f "} - 0(At *) ( 8) 

where all coefficients appear at time n. 


This scheme is now rewritten in the original variables of (1) and (2), 
and centered differences are used to evaluate the divergence terms. Thus, the 
following system of equations is obtained for each interior point of the com- 
putational grid: 


where 


ri +/i At(D,u + DyV) 



*iAt(D,pD, + D,pDy)“ 
1 +pAt(uD, + vD,'_ 






(9) 


A( )-( r‘-( )‘. 
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Approximate Factorization 

The matrix on the left-hand side of (9) can be approximately factored 


as 


'1+MAtDyV 

. l+MAtvD,J 


” 1 + M At D,u D,pD, ~j 

0 1 + p At uD, J 


( 10 ) 


However, this scheme is not unconditionally stable. A small modification of 
(10), which introduces error of order AT*, does result in an unconditionally 
stable scheme. This modification consists of replacing the lower right-hand 
term of the second mutrix by 

1+pAtuD, -^i*At^(a^^p)D,pD,. 

The resulting alternating sweeps of this stable scheme are 


• x-sweep: 


l+pAtD,v pAtD,pDy‘ 
pAt^^^ 1 + pAtvD, 


iriHil 


( 11 ) 


• y- 8 weep: 

l-)-pAtD,u MAtD,pD, 

0 1 +pAt uD,-(p2At’aVp)D,pD, 


|Ap 

|A0 




( 12 ) 


Artificial Viscosity 

To capture shocks, a’^fidal- viscosity terms of the Jameson type^® are 
added to the upper-left portions of each equation and to rl. Respectively, 
these terms are 


- 1 AtAyOj 
-cAIAxDJ 
€AI(AxDJp + A>’DJ p). 

where c is a constant, typically equal to 0.15. 


( 13 ) 
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Alternate achenea for InoorporatlnK artificial viacoaity in a more aelective 
manner are poaaible. One auch alternative ia to awitch the coefficient c from 
a low value in regiona where the flow la aubaonic to a higher value in auper- 
aonic regiona; i.e. , 

c max ^0. , 1 - ^ , (14) 

where c, = artificial viacoaity in ’’aubaonic" zonea 

= maximxim amount of artificial viacoaity added to "auperaonic" 

zonea (M >M^) 
c 

M = local Mach number 

= cutoff Mach number (typically 0.95). 

Another acheme la to vary the viacoaity based on gradients in the flow field; 

e . g. t 


€ -Cq+C, • 


I i+2 1+ I I 



(15) 


When such schemes are used, however, care must be taken to retain conserva- 
tion form in the basic equations. Thus, for example, the first term of Expres- 
sions (13) must be replaced by 


-AtAyD,(cD,). (16) 

In the present work, only the constant (non-switched) form of artificial 
viscosity has been used. Future efforts coiild easily include the switching 
concept . 

Airfoil-Adapted Coordinate Transformation 

To conveniently impose the airfoil boundary conditions, a time-varying 
ooordinute transformation can be introduced to (1) and (2). The general form 
of such a transformation is 

(-((x,y,t), T)"»?(x,y,t), T-t • (17) 
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Equations (1) and (2) tranafonn as 



0T h • 0. 


( 18 ) 


where J is the Jacobian of the transformation and U and V are the contra- 
variant velocity components, given by 

U«e, + 4.u + 4,v 

V-T),+T),U + nyV (19) 

To preserve a factorable form of the equations, cross terms arising when (19) 
is introduced into (18) are explicitly differenced. This technique is illus- 
trated below. 

For the particular sample problem used in this study (see subsection, 
Sample Problem), a simple time-varying sheared-rectilinear coordinate trans- 
formation Is used 


4 -X, Tj-y-S(x,t), T-t, (20) 

where S is the instantaneous airfoil surface location as noted in Pig. 1. For 
this case, the determinant of J is 1 and (19) become 

U» 

V»-<i>,S^ + (i>,(l+s5)-ST. (21) 


Substitution into (18) gives 

Pt + (p<<>()( + (pI -St + U o,Dfl“(ps^o,P^ + (ps,<i),)„ 

pT-l 

1)m 1 ' i 0^) -h.»S,0t0, (22) 
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RtO i7»2-001W 


Figure 1. — Shearing Transformation 


The implicit algorithm developed previously is now applied to this system; 
however the cross terms, grouped on the right-hand side of (22), are handled 
explicitly. The resulting equations are 


U |-s^.7l) AP At(^ 

A<^ + ^AT(aVp)Ap + MAT^^^+V^^ Ad> +ATh = 0. 


Ad +ATh = 0. 


where "a" is the local speed of sound and 

a=0,, V = (l+Sj)d„. )5 = (l + s;)p. (24) 

Replacing the indicated partial derivatives by central-difference operators, 
D, one obtains the following matrix equation: 
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1 ■••M AT(D,n + D, (? - St» AT(D,pD, + D,) 

j 1 + M AT(HD, + VD^ 




Following the factorization procedure previously outlined and Introducing the 
artificial viscosity terms, one obtains the final difference equations: 


* ^ - sweep 

1 + p AT D„(r -Sj)-€ AT Atj Oi| 

■"(7) 


pATD„ pD, 
1 +pATrD„ 




(26) 


where 


I rl |D,pU + D,/)V-t(ACD]+AT?D5)p| 


• ii - sweep 


1 +pATD,n-tATA4Di 


M PD| 


1+pATHD, 






t27) 


Discretization 

In (26) and (27), differences will be centered about (i.j) and will span 

2 mesh widths (3 points); consequently, the operators (D pD ) and (D^pD.) 

n ^ s 

will spread to 4 mesh widths (5 points). The equation for the j-th mesh 
point in (28) may then be written In the form 
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1A){AW},.2^-(B1{AW}h-*-(C1{aw}, +[D){aw},«, +(El{AW}^-{r} 


(28) 


where the coefficients of (AW) for mesh points outside the rang^e j - 2 to j 2 
are sero. If (26) is multiplied by 2(An)j/AT, the coefficients in (28) become 


(A) * 


0 

0 


(P/2Atj)j.j 

0 



(P/2Atj),*i ■ 
0 



(29) 


A similar set of matrices can be obtained from (27). 

Stretchings 

For computational efficiency, simple grid stretchings also are introduced: 

4 = 4(0 and ^ *»?(»)). (^®^ 

The effect of these transformations on the equations to be solved is simply to 
introduce the derivatives of the stretching functions as multipliers of the 
terms containing spatial derivatives. The details of the particular stretchings 
used are discussed in the Sample Problem subsection. 

Solution Procedure 

If equations similar to (28) are written for all mesh points and combined 
into one system of equations corresponding to (26), the result is a system of 
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block-flve-diagonal matrix equations. This can be solved efficiently by a 
splitting procedure. Thus, (26) is rewritten as 

I Ml {aw} -{r}, ( 31 ) 

where 



C, 

D, 

Et 

0 

0 

0 

0 • • 

Bj 

Cs 

Dt 

Es 

0 

0 

0 • • 

As 

Bs 

Cs 

D 3 

E, 

0 

0 • 

0 

• 

A| 

• 

B 4 

• 

C 4 

• 

D 4 

• 

E 4 

• 

0 • • 

• e • 

• 

• 

• 

a 

• 

• 

tea 


or in split form as 

(Ll-(Ul{AW}-{r}, 


where 


”>i 

0 

0 

0 • • " 


” 1 

-6, 


0 0 • • 


>2 

0 

0 • • 


0 

1 

-62 

— «2 0 • 

^>3 

Pi 

rs 

0 • • 


0 

0 

1 

-63 -‘s * • 

0 

“4 

Pa 

>4 • • 

. (ul = 

0 

0 

0 

1 -64 • • 

a 

a 

a 

a a a 


a 

a 

a 

a a a a 

a 

a 

a 

a a a 


a 

a 

a 

a a a a 

. 





_ 



- 


Expanding the product LU and comparing with M. one obtains the following 
relationships between the elements of M , L, and b: 

■yj*Cj+Qj€j_2 + 
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( 33 ) 


(34) 

where the first solution is obtained by a forward sweep in j and the second 

by a reverse sweep. An identical procedure. is used to obtain a solution for 
eons 

V M • / • 

To initialize and terminate the sweeps, the difference equations are mod- 
ified to incorporate the presence of the boundaries. At the lower boundary, 
for example, we choose to backward difference the p- and v-type lorms (so 
that no physical variables are defined outside the computational space) and 
to retain central differences in terms by introducing a i*ow of dummy points 
just below T) = 0. Before and after the sweeps are completed, the boundary 
conditions (discussed below) are enforced. 

Boundary Conditions 

On the airfoil surface, tangential flow is maintained. In the transformed 
coordinates, this is equivalent to 

( 35 ) 

V = 0, £orij = 0 and 

where subscripts LE and T£ denote the airfoil leading and trailing edges, res- 
pectively. For the nonlifting symmetric-airfoil problem studied, the flow is 
symmetric about n= 0; consequently, only half the flow field is modeled. 

Thus , 

forT )«^0 and ( ^ or ( ^ • (36) 

The uniform onset flow in the far field is assumed to be undistui'bed : i.e. , 
acoustic waves originating at the airfoil do not have sufficient time to reach 
the outer boundaries. Thus, 

(u,v)-=U,0). 


Then, the solution is obtained by sequentially solving 

|U|{AW}-{2f. 

Thus, 

AW, - z, + 6,(AW,„) +€ ,(AW,*,), 


13 



or 

v« <^ 1 , -0, 

u ■ - S, ^ ■ 1 . 

Hence . 


forTj«Tj^ and (37) 


where the subscripts "max" denote the far-field boundaries of the compu- 
tational mesh. In Fig. 2, the boundary conditions are marked on a sketch of 
the computational space. 



Figure Z — Boundary Conditiont in the Computational Grid. 


Implementation 
(Codes UFL03 and UFL04) 

The algorithm was coded for the pulsating-airfoil sample problem to be 
described below. A version (UFL03) using the original untransformed equa- 
tions was written, as well as a version (UFL04) employing the time-varying 
coordinates. For UFL03, the airfoil boundary conditions were applied on the 
slit, y = 0 (mean chordline); whereas, in U FLO 4, they were applied on the 
r, = 0 line coincident with the instantaneous airfoil surface position. 
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5 ' RESULTS 


Sample Problem 

A problem that has been analysed by several researchers is that of com- 
puting unsteady pressures on a parabolic-arc airfoil which successively 
thickens and thins during its travel, as shown in Fig. 3. The Mach number 
for this example is 0.85. The equations governing the variation of thickness 
are 

TH»0, 30iT, 

where TH is the midcord thickness ratio and T is time, measured in chord- 
lengths traveled. Consequently, the airfoil initially has zero thickness, grows 
to its maximum thickness of 10% after traveling 15 chordlengths and returns to 


S Ik. tl - THIt) • SgU) 

Sglx) 10% THICK BICIRCULAR ARC AIRFOIL 



L CHORDLENGTHS TRAVELED 

nto-irsaoojw 


Figure 3. — Puluting Parabolic-Arc Airfoil. 
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lero thickness after traveling a toUil of 30 chordlengths. During the course of 
this travel, the variation of thickness causes an interesting flow-structure. 

A strong shocV wave forms on the airfoil ai it thickens; subsequently, as the 
airfoil thins, the shock propagates rapidly upstream and leaves the airfoil 
nose to enter the oncoming flow. The numerical computation of this extensive 
shock motion is a rigorous test for unsteady transonic aerodynamic codes. It 
might be noted that, because of their basic theoretical limitations, the methods 
of Ref. 1-5 are unable to handle such cases of large shock motion. 

The computational grid used consists of 152 points in the streamwise di- 
rection and 40 in the stream-normal direction. To facilitate compariecns, the 
grid is patterned after that of Ref. 15. In the streamwise direction, the 
grid is uniform over the Interval that extends from one chordlength upstream 
of the airfoil nose to the trailing edge: to either side of this interval, the grid 
is smoothly stretched to the boundaries located more than 30 chordlengths from 
the airfoil. In the stream-normal direction, the grid is uniform from the air- 
foil surface to a distance of 0.2 chordlengths; beyond this point, the grid is 
stretched smoothly to a boundary also more than 30 chordlengths from the air- 
foil. The minimum grid spacing is roughly 0.02 chordlengths in each direc- 
tion. From studies of grid variation, it was concluded that the solution is 
sencitivA to the choice of grids but that the present choice is adequate because 
it combines a fine- grid structure near the airfoil with boundaries sufficiently 
far removed for the present calculations. Appendix A presents the actual 
mesh used. 

Unless otherwise noted, all calculations using either UFL03 or UFL04 were 
performed with a time step of AT = 0.02 chordlengths traveled and an artifical 
viscosity coefficient of c = 0.15. Runs varying these parameters are described 
in separate subsections below. 

UFL0 3 and UFL04 Results 

Using UFL03, a time history of the flow about the pulsating airfoil was 
calculated. The resulting pressure coefficient distributions for six time slices 
are shown in Fig. 4. At the first time slice, the flow is subcritical. During 
the next two, a shock forms, strengthens, and moves aft. A slight re-expan- 
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■ion occurring behind the shock can be seen in Fig. 4, Sheet 1. It should be 
noted that a significant lag occurs between the time that the airfoil reaches 
maximum thickness (T = 15) and the point at which maximum shock strength 
is attained (T = 18.25). In the next three time slices (note that different 
scales are used), the shock moves rapidly forward, while diminishing in 
strength, and leaves the airfoil. 

To determine the effect of applying the boundary conditions on the slit 
rather than the actual airfoil surface, results of UFL03 and UFL04 are com- 
pared in Fig. 5. To obtain a more dramatic difference, a 15% airfoil was used 
in place of the 10% airfoil previously studied; consequently, different time 
slices are shown. At the early time shoes (Fig. 5, Sheet 1 and 2), during 
which the shock is formed, the results are practically identical. At later times 
noticeable differences occur. Compaf*;ng time T = 15 (Fig. 5, Sheet 2) with 
T = 25 (Fig. 5, Sheet 3) and time T = 25 with T = 30 (Fig. 5. Sheet 3), one 
sees that the shock speeds computed by UFLO 3 during these time intervals are 
greater than those computed by UFL04 (in which the boundary conditions ar^' 
correctly apphed on the airfoil surface). A comparison of time T = 30 with T = 
35 indicates thet this trend persists even after the airfoil has returned to 
zero thickness. 


Comparison With Other Methods 


Only one other operational code (Ref. 15) exists that solves the un- 
steady-potential-flow equations in conservation form with boundary conditions 
correctly applied on the instantaneous airfoil surface. (The method of Ref. 14 
uses a primitive-variable formulation rather than introducing the potential 
function itself. Consequently, vorticity is numerically created in the flow 
field and degenerates the solution accuracy, particularly behind strong 
shocks.) Thus, comparisons between results of the present method and Ref. 15 
serve as a check on the correctness of both. Furthermore, since Ref. 14 and 
15 give comparisons with other, less exact formulations (Ref. 10 and 13), 
we will omit such investigations here. 
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UFL03& UFL04 
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Fi^rt 5. - Comparition of UFL03 and UFL04 foe a 15%-Thick Airfoil. (Shaat 1 of 3). 


UFL03& UFL04 
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Figur* 5. - Comparison of UFL03 and UFL04 for a 15%- Thick Airfoil, (Shaat 2 of 3), 





For the 10%-thick puldatin^ airfoil. res\ilts of UFL04 and Goorjian's 
method (Ref. 15) are shown in Fig. 8. The UFL04 results are so labeled In the 
plot legend; the Goorjian results are labeled "Ref. 15". As can be seen, 
agreement in generally good; however, 1. _re are several minor discrepancies. 
During shock formation, UFL04 predicts a slightly weaker, farther aft shock. 
Since the shock- location difference is only one mesh width (see Fig. 6, Sheet 
3). it is considered insignificant. During shock propagation, the Goorjian- 
method shock is smeared over larger distances and attenuates more rapidly 
than that of UFL04. The better shock resolution of the present method is 
considered to be due to the use of artifical viscosity in lieu of artificial com- 
pressibility as used by Goorjian. Some of the remaining disparity is due to 
the use of slightly different grids in the two methods.* Also, the calculations 
were performed with different time steps. The effect of this parameter in 
UFL04 is discussed in a following subsection in which it is established that 
the time step used is adequate for accurate solutions; consequently, only 
small discrepancies can be attributed to this factor. In summary, the com- 
parisons show that the present method and that of Ref. 15 give almost equiv- 
alent results for the sample problem. The only significant difference in the 
solutions is the degree of shock smearing and, in this aspect, the present 
program appears to perform somewhat better. 

Stability, Accuracy, and Computing Time 

To determine the effect of the time step on the stability and accuracy 
of the difference scheme, both UFL03 and UFL04 were rerun varying this 
parameter. The results of these investigations are summarized in Fig. 7 and 
8(UFL03) and Fig. 9 and 10 njFL04) , where the computed pressures are 
compared at two time slices. InUFL04, the time step was successfully varied 
from 0.01 to 0.20 chordlenghts traveled without introducing numerical insta- 
bility. Larger time steps were not attempted. In UFL03,the time step was 
varied from 0.01 to 0.065. Since trends similar to those of UFL04 were 
found, the larger step, 0.20, was ;,ot run. It is concluded that the algorithm 
has excellent stability and that the introduction of the time- varying coordinates 
(UFL04) does not degrade this quality. 
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From inspection of Fig. 7 through 10, conclusions on the effect of time 
step on solution accuracy can be made . As can be seen . reasonable results 
were obtained for AT as large as 0.033. With AT = 0.065, however, the com- 
puted shock speed during propagation is significantly less than it should be 
and, with aT = 0.20, the results correlate poorly altogether with those of 
the other runs. 

As mentioned previously, the results of Fig. 6 were obtained with aT = 
0.020. Figures 16 and 18 show that the effects of using a smaller step are 
quite small. Thus, the step used in the comparison studies was adequate. 

Studies discussed in the next subsection indicate that a smaller amount 
of artificial viscosity can be successfully used. Since, for non-rero mesh 
widths, the artificial viscosity terms in (26) and (27) can be thought of as 
error terms of order AT, any reduction in these terms should improve accu- 
racy for a given step size. Thus, it is likely that, by reducing the amount 
of artificial viscosity, accurate solutions could be obtained with a somewhat 
larger step size than the 0,033 discussed above. 

The computational time required for the algorithm is roughly 7.0 x 10~5 
seconds per time step per grid point on the CDC 7600 computer. For the sam- 
ple problem with aT = 0.033, the time required for a time history of 32 chord- 
lengths tra' eled is about 6 minutes. 

Effect of Artificial Viscosity 

To determine the amount of artificial viscosity necessary to sharply re- 
solve shocks, variations of the parameter e were made in the sample problem. 
The results of this investigation are shown in Fig. 11 for a single point in the 
time history. It can be seen from this figure that the value 0.15 used in all 
previous calculations can be vairied by approximuteiy 20% without introducing 
exce.ssive overshoots or smearing near the shock. As expected, variations of 
c have no effect away from the shock. 
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Figure 6. — Comp8ri:on of UFL04 erith Goorjian Program (Sheet 1 of P). 
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Figure 6. - Comparison of UFL04 with Goorjian Program (Sheet 2 of 6) 
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Figure 6. — Comparison of UFL04 with Goorjian Profpvn (Shaat 3 of 6). 
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Figure 6. — Comparison of UFL04 with Goorjisn Program (sheet 6 of 6) 
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Figure 8. - Fffect of Integration Step: UFL03 at T » 3Z0. 
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6 - CONCLUSIONS AND RECOMMENDATIONS 


An accurate algcriihm with excellent numerical stability has been devel- 
oped for the two-dimensional unsteady full-potential equation in conservation 
form, using time-varying sheared-rectilinear , body-conforming coordinates. 
The method has been demonstrated on the highly nonlinear transonic flow 
arising from an airfoil with pulsating thickness. The results of this demon- 
stration have been verified by comparison with those reported in Ref. 15. 

Two recommendations are made to increase the accuracy and efficiency 
of the procedure. The first is to introduce switched artificial viscosity as 
discussed in the THEORY section. This change would localize the artificial- 
viscosity terms to regions where they are needed, thereby improving the ac- 
curacy for a given time step. Consequently, larger time steps could be used 
with a resultant decrease in computing time. The second suggestion is to 
modify the manner in which the operators D^pD, and in (26) and 

I 27) are evaluated. A three-point scheme could be constructed in place of a 
five-point scheme. To accomplish this, densities midway between mesh points 
would be evaluated by averaging the values from the two adjacent nodes. 

This modification would result in reducing the system in (31) from block-fivc- 
diagonal to block-three-diagonal. It can be shown that the computing time 
required to solve the smaller-bandwidth system is roughly half that required 
for the larger system. 

To enable the method to handle realistic blunt-nosed airfoils, it is further 
recommended that the code be modified to use a curvilinear coor'^inate trans- 
formation. In particular, a parabolic (C-mesh) transformation is suggested. 
Additional efforts should then be focused on extensions to the case of a lift- 
ing airfoil and improv’ed treatment of the far-field boundary conditions. 
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